D-place FARM documentation: Module 3

library(png)

format(Sys.time(), format="%d_%b_%Y")

[1] "27_Feb_2017"

## First consolidate the available files into a single table
    
      path <- "~/Box Sync/Four model compare/Module 2 extinct"
           
     
           setwd(path)

Error in setwd(path) : cannot change working directory

load('~/Box Sync/colliding ranges/Simulations_humans/Results/available daily summaries/Four_model_compare_results_02_Mar_2017_crop_to_3481.Rdata')
extant <- Concatenated_data
extant

setwd("~/Box Sync/colliding ranges/Simulations_humans/Results/available daily summaries")

The working directory was changed to /Users/Ty/Box Sync/colliding ranges/Simulations_humans/Results/available daily summaries inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the the working directory for notebook chunks.

details <- file.info(list.files())
trimmed_details <- details[which(list.files() == list.files(pattern = "Four_model_compare_results_extinct*")),]
ord <- order(trimmed_details$mtime, decreasing = TRUE)
rownames(trimmed_details[ord,])[1]

[1] "Four_model_compare_results_extinct_10_Mar_2017.Rdata"

load(rownames(trimmed_details[ord,])[1])
extinct <- Concatenated_data
trimmed_details <- details[which(list.files() != list.files(pattern = "Four_model_compare_results_extinct*")),]
ord <- order(trimmed_details$mtime, decreasing = TRUE)
rownames(trimmed_details[ord,])[1]

[1] "Four_model_compare_results_10_Mar_2017_crop_to_4897.Rdata"

load(rownames(trimmed_details[ord,])[1])
extant <- Concatenated_data

dim(extinct)

[1] 25492    62

dim(extant)

[1] 14660    62

for(i in c(9,10,11,12,14,15,16,17,19,20,21,22,24,25,26,27,29,30,31,32)){
    extinct[which(is.nan(as.numeric(as.character(extinct[, i]))) == TRUE), i] <- NA
}
for(i in c(9,10,11,12,14,15,16,17,19,20,21,22,24,25,26,27,29,30,31,32)){
    extant[which(is.nan(as.numeric(as.character(extant[, i]))) == TRUE), i] <- NA
}
i <- 19
for(i in c(20,21,24,25,26,27)){
    extinct[which(as.numeric(as.character(extinct[, i])) == 0), i] <- NA
}
for(i in c(20,21,24,25,26,27)){
    extant[which(as.numeric(as.character(extant[, i])) == 0), i] <- NA
}
xlimit <- c(0,1)
ylimit <- c(0,600)
maincex <- 0.9
png(file="Global_success_rate_per_parameter.png", width=8.5, height=11, units="in", res=300)
par(mfrow=c(5,4), mar=c(3,3,3,0))
hist(as.numeric(as.character(extinct[,9])), main="speciation of F in F env", col=adjustcolor("firebrick", alpha=0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,9])), main="speciation of F in F env", col=adjustcolor("cornflowerblue", alpha=0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,10])), main="speciation of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[,10])), main="speciation of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,11])), main="speciation of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,11])), main="speciation of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[,12])), main="speciation of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[,12])), main="speciation of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
#######
hist(as.numeric(as.character(extinct[, 14])), main="extinction of F in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 14])), main="extinction of F in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 15])), main="extinction of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 15])), main="extinction of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 16])), main="extinction of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 16])), main="extinction of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 17])), main="extinction of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 17])), main="extinction of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
######
hist(as.numeric(as.character(extinct[, 29])), main="arisal of F in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 29])), main="arisal of F in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 30])), main="arisal of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 30])), main="arisal of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 31])), main="arisal of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 31])), main="arisal of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 32])), main="arisal of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 32])), main="arisal of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
######
hist(as.numeric(as.character(extinct[, 19])), main="NOPE -- Diffusion: source F, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex)
hist(as.numeric(as.character(extant[, 19])), main="NOPE -- Diffusion: source F, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 20])), main="Diffusion: source D, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 20])), main="Diffusion: source D, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 21])), main="Diffusion: source F, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 21])), main="Diffusion: source F, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 22])), main="NOPE -- Diffusion: source D, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex)

hist(as.numeric(as.character(extant[, 22])), main="NOPE -- Diffusion: source D, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex, add=TRUE)
####
hist(as.numeric(as.character(extinct[, 24])), main="Takeover: source F, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 24])), main="Takeover: source F, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 25])), main="Takeover: source D, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 25])), main="Takeover: source D, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 26])), main="Takeover: source F, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 26])), main="Takeover: source F, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
hist(as.numeric(as.character(extinct[, 27])), main="Takeover: source D, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)

hist(as.numeric(as.character(extant[, 27])), main="Takeover: source D, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)
dev.off()

null device 
          1 

png(file="extiction minus extant per outcome.png", width=8.5, height=11, units="in", res=300)
par(mfrow=c(3,1))
plot(as.numeric(as.character(extinct[,9])), as.numeric(as.character(extinct[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("firebrick", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))
plot(as.numeric(as.character(extant[,9])), as.numeric(as.character(extant[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("cornflowerblue", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))

plot(as.numeric(as.character(extinct[,9])), as.numeric(as.character(extinct[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("firebrick", alpha=0.2), pch=19, cex=0.6, ylim=c(0,1))
points(as.numeric(as.character(extant[,9])), as.numeric(as.character(extant[,14])), xlab="speciation", ylab="extinction", col= adjustcolor("cornflowerblue", alpha=0.2), pch=19, cex=0.6)
dev.off()

null device 
          1 

xfit1 <- seq(xmin, xmax, length= 100)

Error in is.finite(from) : 
  default method not implemented for type 'closure'

load('~/Box Sync/colliding ranges/Simulations_humans/Available trees/real.analysis.RData')
Concatenated_data <- extant
#Concatenated_data <- Concatenated_data[Concatenated_data[, 2] == "stats.no.bTO", ]
#Concatenated_data <- Concatenated_data[Concatenated_data[, 6] != "05", ]
# Concatenated_data[, 6] <- as.numeric(Concatenated_data[, 6])
# # Concatenated_data[original[, 2] == "background_takeover", 6] <-  Concatenated_data[original[, 2] == "background_takeover", 6] + 4
Concatenated_data[, 6] <- factor(Concatenated_data[, 6])
#head(Concatenated_data)
#names(Concatenated_data)
PCAdata <- Concatenated_data[, -(1:35)]
PCAdata <- PCAdata[, -12]
PCAdata <- apply(PCAdata, 2, as.numeric)
remove <- apply(is.na(PCAdata), 1, any)
PCAdata <- PCAdata[!remove, ]
# Predictions
library(randomForest)

randomForest 4.6-12
Type rfNews() to see new features/changes/bug fixes.

data.analysis.comp2 <- data.frame("Model" = as.factor(Concatenated_data[!remove, 6]),
                                  PCAdata)
data.analysis.comp2$sprate <- data.analysis.comp2$trait_1_speciation/data.analysis.comp2$trait_2_speciation
data.analysis.comp2$extrate <- data.analysis.comp2$trait_1_extinction/data.analysis.comp2$trait_2_extinction
#load("Real_phy/real.analysis.RData")
a <- as.data.frame(real.analysis$results_summary_of_single_value_outputs)
a$sprate <- a$trait_1_speciation / a$trait_2_speciation
a$extrate <- a$trait_1_extinction / a$trait_2_extinction
data.analysis.comp3 <- data.analysis.comp2[, -c(2, 13:14, 16:20)]
#data.analysis.comp3 <- data.analysis.comp3[data.analysis.comp3$Model %in% 1:4, ]
#data.analysis.comp3$Model <- factor(data.analysis.comp3$Model)
#sub <- unlist(lapply(as.list(c(1:4)), function(x, y) {
#  sample(which(y$Model == x), min(table(data.analysis.comp3$Model)))},
#  y = data.analysis.comp3))
# data.analysis.comp3 <- data.analysis.comp3[sub, ]
fun <- function(x, y, per = .33) {sample(which(y$Model == x), round(table(y$Model)[1]*per))}
sub.test <- unlist(lapply(as.list(paste0(0, c(1:4))), fun,
                          y = data.analysis.comp3))
test2 <- data.analysis.comp3[sub.test, 2:ncol(data.analysis.comp3)]
test1 <- data.analysis.comp3[sub.test, 1]
train <- data.analysis.comp3[-sub.test, ]
for(i in 1:100){
(fit <- randomForest(Model ~ ., data=train, xtest = test2, ytest = test1, 
                    importance=TRUE, ntree=1000, keep.forest = TRUE, replace=FALSE))
predictions <- predict(fit, 
                       a,
                       type="prob")
predictions
save(fit, file=paste0("~/Box Sync/colliding ranges/Simulations_humans/Results/RF_daily_output/RF_daily_output_", format(Sys.time(), format="%d_%b_%Y"), "_",i, "_NoREPLACEMENT_.Rdata"))
}

save(fit, file=paste0("~/Box Sync/colliding ranges/Simulations_humans/Results/RF_daily_output/RF_daily_output_", format(Sys.time(), format="%d_%b_%Y"),".Rdata"))

plot(fit, ylim=c(0,1))

labs <- c("Basic", "+Diffusion", "+Takeover", "+Diffusion +Takeover")

# bar plot
png("Prob_aus.png", width = 25, height = 25, res = 300, units = "in")
par(mar = c(8, 8, 1, 1))
pred <- setNames(as.numeric(predictions), labs)
cols <- rev(c("darkgreen", "red", "blue", "darkorange1"))
barplot(pred, col = cols, ylab = "Proability", cex.lab = 3, cex.names = 2)
dev.off()

null device 
          1 

# Plot confusion matrix
png("Conffusion_matrix_all.png", width = 25, height = 25, res=300, units="in")
par(mar = c(10, 11, 1, 1))
colors1 <- colorRampPalette(colors = c("#f0f0f0", "#bdbdbd","#636363"))
prop <- apply(fit$confusion[, -5], 2, function(x){x / sum(x)}) * 100
image(prop, col = colors1(20), axes=FALSE)
axis(1, at=c(0, .33, .66, 1), labels=labs, tick = FALSE, line = FALSE, cex.axis = 3.5, pos = -.19)
axis(2, at=c(0, .33, .66, 1), labels=labs, tick = FALSE, line = FALSE, cex.axis = 3.5)
mtext("ACTUAL", side = 1, padj = 3, cex = 4)
mtext("PREDICTED", side = 2, padj = -3, cex = 4)
for(i in 1:4) {
  for(j in 1:4) {
    text(x = c(0, .33, .66, 1)[i], y = c(0, .33, .66, 1)[j], paste0(round(prop[i, j], 2), "%"),
         cex = 5)
  }
}
dev.off()

null device 
          1 

importance(fit)

                                                   01        02        03
Pylo_diversity_is_sum_of_BL                  36.68603 25.692246  44.88765
average_phylogenetic_diversity_is_mean_of_BL 43.12789 34.908001  47.98480
variance_Pylo_diversity_is_variance_of_BL    38.00313 20.262707  19.18043
F_quadratic_entropy_is_sum_of_PD             28.49584 23.877597  40.27946
Mean_pairwise_distance                       35.73981 27.249223  32.09020
variance_pairwise_distance                   39.27909 38.069218  40.92584
Evolutionary_distinctiveness_sum             37.35990 28.081845  45.31515
mean_Phylogenetic_isolation                  41.20750 35.209163  48.94046
variance_Phylogenetic_isolation              35.88075  6.978819  28.67514
gamma                                        27.59011 28.762859  76.57843
extinction_per_speciation                    11.04195 51.084694  65.78277
transition_from_trait_1_to_2                 28.53323 27.846818  47.04102
transition_from_trait_2_to_1                 36.78696 30.942883  50.69606
transition_rate_ratio_1to2_over_2to1         54.87487 38.192677  38.58033
Phylogenetic_signal                          66.29483 45.359260  70.57335
spatial.tests.fora                           53.03264 95.790134  98.20864
spatial.tests.dom                            86.83866 40.905940  53.51233
prevalence                                   55.09973 15.677005  39.96226
sprate                                       57.87711 32.873304 105.88921
extrate                                      29.72153 29.191046  40.11454
                                                   04 MeanDecreaseAccuracy
Pylo_diversity_is_sum_of_BL                  45.89839             67.10212
average_phylogenetic_diversity_is_mean_of_BL 35.72290             71.53199
variance_Pylo_diversity_is_variance_of_BL    31.31532             59.79845
F_quadratic_entropy_is_sum_of_PD             42.73530             76.03108
Mean_pairwise_distance                       42.95949             69.53536
variance_pairwise_distance                   37.44388             76.20714
Evolutionary_distinctiveness_sum             46.70827             67.45126
mean_Phylogenetic_isolation                  36.35073             72.84609
variance_Phylogenetic_isolation              23.59933             50.96412
gamma                                        51.86023            104.53469
extinction_per_speciation                    53.99058             87.53765
transition_from_trait_1_to_2                 36.98221             64.93037
transition_from_trait_2_to_1                 42.85645             68.22442
transition_rate_ratio_1to2_over_2to1         51.09836             93.31104
Phylogenetic_signal                          46.37239             91.08957
spatial.tests.fora                           76.18005            138.82996
spatial.tests.dom                            64.49868            119.35092
prevalence                                   64.89918             86.99762
sprate                                       41.98590            109.24396
extrate                                      37.94293             69.14048
                                             MeanDecreaseGini
Pylo_diversity_is_sum_of_BL                          271.5521
average_phylogenetic_diversity_is_mean_of_BL         241.4081
variance_Pylo_diversity_is_variance_of_BL            222.2757
F_quadratic_entropy_is_sum_of_PD                     247.9364
Mean_pairwise_distance                               272.9117
variance_pairwise_distance                           259.1696
Evolutionary_distinctiveness_sum                     277.8165
mean_Phylogenetic_isolation                          241.4813
variance_Phylogenetic_isolation                      220.4116
gamma                                                283.5771
extinction_per_speciation                            340.7512
transition_from_trait_1_to_2                         335.8418
transition_from_trait_2_to_1                         356.5606
transition_rate_ratio_1to2_over_2to1                 288.0278
Phylogenetic_signal                                  418.6853
spatial.tests.fora                                   536.9193
spatial.tests.dom                                    344.7908
prevalence                                           303.3840
sprate                                               342.5898
extrate                                              260.9412

# Variables importance
imp <- importance(fit)
imp <- apply(imp, 2, function(x) (x - min(x))/(max(x) - min(x)))
imp <- imp[sort(imp[, 5], index.return = TRUE, decreasing = TRUE)$ix, ]
names <- rownames(imp)
names[names == "spatial.tests.fora"] <- "Space F"
names[names == "spatial.tests.dom"] <- "Space D"
names[names == "sprate"] <- "Sp(ratio)"
names[names == "transition_from_trait_1_to_2"] <- "TR(1-2)"
names[names == "transition_from_trait_2_to_1"] <- "TR(2-1)"
names[names == "Phylogenetic_signal"] <- "PhySig(D)"
names[names == "Evolutionary_distinctiveness_sum"] <- "EDsum"
names[names == "Pylo_diversity_is_sum_of_BL"] <- "PDsum"
names[names == "transition_rate_ratio_1to2_over_2to1"] <- "TR(ratio)"
names[names == "gamma"] <- "Gamma"
names[names == "mean_Phylogenetic_isolation"] <- "MPI"
names[names == "extrate"] <- "Ext(ratio)"
names[names == "average_phylogenetic_diversity_is_mean_of_BL"] <- "PDmean"
names[names == "extinction_per_speciation"] <- "DR"
names[names == "variance_Phylogenetic_isolation"] <- "VPI"
names[names == "F_quadratic_entropy_is_sum_of_PD"] <- "F"
names[names == "Mean_pairwise_distance"] <- "MPD"
names[names == "variance_Pylo_diversity_is_variance_of_BL"] <- "PDvar"
names[names == "variance_pairwise_distance"] <- "VPD"
png("var_import_all.png", width = 25, height = 25, unit="in", res=300)
par(mar = c(10, 18, 1, 1))
plot(x = rev(imp[, 5]), y = 1:nrow(imp), type = "l", yaxt = "n", 
     ylab = "", xlab = "Variable Importance",
     xlim = c(0, 1), lwd = 2, cex.lab = 4)
for (i in 1:nrow(imp)) {
  abline(h = i, lty = 3, col = "gray80")
}
abline(v = seq(0, 1, 1/19), lty = 3, col = "gray80")
lines(x = rev(imp[, 4]), y = 1:nrow(imp), col = "darkgreen", lwd = 2)
lines(x = rev(imp[, 3]), y = 1:nrow(imp), col = "red", lwd = 2)
lines(x = rev(imp[, 2]), y = 1:nrow(imp), col = "blue", lwd = 2)
lines(x = rev(imp[, 1]), y = 1:nrow(imp), col = "darkorange1", lwd = 2)
lines(x = rev(imp[, 5]), y = 1:nrow(imp), lwd = 3)
points(x = rev(imp[, 4]), y = 1:nrow(imp), col = "darkgreen", cex = 2)
points(x = rev(imp[, 3]), y = 1:nrow(imp), col = "red", cex = 2)
points(x = rev(imp[, 2]), y = 1:nrow(imp), col = "blue", cex = 2)
points(x = rev(imp[, 1]), y = 1:nrow(imp), col = "darkorange1", cex = 2)
points(x = rev(imp[, 5]), y = 1:nrow(imp), pch = 20, cex = 3)
text(y = 1:nrow(imp), x = par("usr")[1] - .17, labels = rev(names),
     srt = 0, pos = 4, xpd = T, cex = 4)
dev.off()

null device 
          1 

par(mfrow=c(2,3))
# Box plots
boxplot(spatial.tests.fora ~ Model, data = data.analysis.comp3)
abline(h = a$spatial.tests.fora, col = "red", lty = 2)

boxplot(spatial.tests.dom ~ Model, data = data.analysis.comp3)
abline(h = a$spatial.tests.fora, col = "red", lty = 2)

boxplot(log(sprate) ~ Model, data = data.analysis.comp3, ylim = c(-10, 10))
abline(h = log(a$sprate), col = "red", lty = 2)

boxplot(log(extrate) ~ Model, data = data.analysis.comp3, ylim = c(-10, 10))
abline(h = log(a$extrate), col = "red", lty = 2)

boxplot(log(transition_rate_ratio_1to2_over_2to1) ~ Model, data = data.analysis.comp3)
abline(h = log(a$sprate), col = "red", lty = 2)

boxplot(Phylogenetic_signal ~ Model, data = data.analysis.comp3, ylim = c(0, 1))
abline(h = a$Phylogenetic_signal, col = "red", lty = 2)

str(fit)

List of 19
 $ call           : language randomForest(formula = Model ~ ., data = train, xtest = test2, ytest = test1,      importance = TRUE, ntree = 2400, keep.forest = TRUE)
 $ type           : chr "classification"
 $ predicted      : Factor w/ 4 levels "01","02","03",..: 1 3 1 1 3 3 3 1 3 1 ...
  ..- attr(*, "names")= chr [1:8885] "1" "2" "3" "4" ...
 $ err.rate       : num [1:2400, 1:5] 0.547 0.547 0.549 0.535 0.533 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:5] "OOB" "01" "02" "03" ...
 $ confusion      : num [1:4, 1:5] 1392 298 434 90 230 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "01" "02" "03" "04"
  .. ..$ : chr [1:5] "01" "02" "03" "04" ...
 $ votes          : matrix [1:8885, 1:4] 0.549 0.346 0.906 0.902 0.22 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:8885] "1" "2" "3" "4" ...
  .. ..$ : chr [1:4] "01" "02" "03" "04"
  ..- attr(*, "class")= chr [1:2] "matrix" "votes"
 $ oob.times      : num [1:8885] 907 876 855 911 850 863 879 836 883 902 ...
 $ classes        : chr [1:4] "01" "02" "03" "04"
 $ importance     : num [1:20, 1:6] 0.0248 0.0262 0.0149 0.012 0.0186 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  .. ..$ : chr [1:6] "01" "02" "03" "04" ...
 $ importanceSD   : num [1:20, 1:5] 0.000668 0.000584 0.00033 0.000355 0.000494 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  .. ..$ : chr [1:5] "01" "02" "03" "04" ...
 $ localImportance: NULL
 $ proximity      : NULL
 $ ntree          : num 2400
 $ mtry           : num 4
 $ forest         :List of 14
  ..$ ndbigtree : int [1:2400] 3853 3761 3829 3805 3823 3977 3887 3917 3887 3907 ...
  ..$ nodestatus: int [1:4035, 1:2400] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ bestvar   : int [1:4035, 1:2400] 13 8 11 16 7 17 10 7 18 11 ...
  ..$ treemap   : int [1:4035, 1:2, 1:2400] 2 4 6 8 10 12 14 16 18 20 ...
  ..$ nodepred  : int [1:4035, 1:2400] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ xbestsplit: num [1:4035, 1:2400] 193.4137 0.0287 0.9995 25.8679 33.4238 ...
  ..$ pid       : num [1:4] 1 1 1 1
  ..$ cutoff    : num [1:4] 0.25 0.25 0.25 0.25
  ..$ ncat      : Named int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "names")= chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  ..$ maxcat    : int 1
  ..$ nrnodes   : int 4035
  ..$ ntree     : num 2400
  ..$ nclass    : int 4
  ..$ xlevels   :List of 20
  .. ..$ Pylo_diversity_is_sum_of_BL                 : num 0
  .. ..$ average_phylogenetic_diversity_is_mean_of_BL: num 0
  .. ..$ variance_Pylo_diversity_is_variance_of_BL   : num 0
  .. ..$ F_quadratic_entropy_is_sum_of_PD            : num 0
  .. ..$ Mean_pairwise_distance                      : num 0
  .. ..$ variance_pairwise_distance                  : num 0
  .. ..$ Evolutionary_distinctiveness_sum            : num 0
  .. ..$ mean_Phylogenetic_isolation                 : num 0
  .. ..$ variance_Phylogenetic_isolation             : num 0
  .. ..$ gamma                                       : num 0
  .. ..$ extinction_per_speciation                   : num 0
  .. ..$ transition_from_trait_1_to_2                : num 0
  .. ..$ transition_from_trait_2_to_1                : num 0
  .. ..$ transition_rate_ratio_1to2_over_2to1        : num 0
  .. ..$ Phylogenetic_signal                         : num 0
  .. ..$ spatial.tests.fora                          : num 0
  .. ..$ spatial.tests.dom                           : num 0
  .. ..$ prevalence                                  : num 0
  .. ..$ sprate                                      : num 0
  .. ..$ extrate                                     : num 0
 $ y              : Factor w/ 4 levels "01","02","03",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "names")= chr [1:8885] "1" "2" "3" "4" ...
 $ test           :List of 5
  ..$ predicted: Factor w/ 4 levels "01","02","03",..: 1 3 1 1 1 3 1 3 2 1 ...
  .. ..- attr(*, "names")= chr [1:188] "1410" "86" "2307" "1834" ...
  ..$ err.rate : num [1:2400, 1:5] 0.559 0.537 0.473 0.479 0.473 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:5] "Test" "01" "02" "03" ...
  ..$ confusion: num [1:4, 1:5] 27 6 10 4 6 24 1 8 11 2 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:4] "01" "02" "03" "04"
  .. .. ..$ : chr [1:5] "01" "02" "03" "04" ...
  ..$ votes    : matrix [1:188, 1:4] 0.94 0.354 0.532 0.543 0.525 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:188] "1410" "86" "2307" "1834" ...
  .. .. ..$ : chr [1:4] "01" "02" "03" "04"
  .. ..- attr(*, "class")= chr [1:2] "matrix" "votes"
  ..$ proximity: NULL
 $ inbag          : NULL
 $ terms          :Classes 'terms', 'formula'  language Model ~ Pylo_diversity_is_sum_of_BL + average_phylogenetic_diversity_is_mean_of_BL +      variance_Pylo_diversity_is_variance_of_BL + F_quadratic_entropy_is_sum_of_PD +  ...
  .. ..- attr(*, "variables")= language list(Model, Pylo_diversity_is_sum_of_BL, average_phylogenetic_diversity_is_mean_of_BL,      variance_Pylo_diversity_is_variance_of_BL, F_quadratic_entropy_is_sum_of_PD,  ...
  .. ..- attr(*, "factors")= int [1:21, 1:20] 0 1 0 0 0 0 0 0 0 0 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:21] "Model" "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" ...
  .. .. .. ..$ : chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  .. ..- attr(*, "term.labels")= chr [1:20] "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" "F_quadratic_entropy_is_sum_of_PD" ...
  .. ..- attr(*, "order")= int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "intercept")= num 0
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(Model, Pylo_diversity_is_sum_of_BL, average_phylogenetic_diversity_is_mean_of_BL,      variance_Pylo_diversity_is_variance_of_BL, F_quadratic_entropy_is_sum_of_PD,  ...
  .. ..- attr(*, "dataClasses")= Named chr [1:21] "factor" "numeric" "numeric" "numeric" ...
  .. .. ..- attr(*, "names")= chr [1:21] "Model" "Pylo_diversity_is_sum_of_BL" "average_phylogenetic_diversity_is_mean_of_BL" "variance_Pylo_diversity_is_variance_of_BL" ...
 - attr(*, "class")= chr [1:2] "randomForest.formula" "randomForest"